#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <string.h>
#include <stdarg.h>
#include <ctype.h>

float hypot3(float x, float y, float z) {return sqrt(x*x+y*y+z*z);}

float dot(float x[3], float y[3]) {return x[0]*y[0]+x[1]*y[1]+x[2]*y[2];}

void *malloc_check(size_t size)
{
  void *p = malloc(size);
  if (p==NULL) {
    fprintf(stderr, "Error allocating %ld bytes\n", size);
    perror("malloc"); exit(EXIT_FAILURE);
  }
  return p;
}

void *calloc_check(size_t size, size_t n)
{
  void *p = calloc(size,n);
  if (p==NULL) {
    fprintf(stderr, "Error allocating %ld * %ld bytes\n", n, size);
    perror("calloc"); exit(EXIT_FAILURE);
  }
  return p;
}

/* allocate a float matrix with subscript range m[0..nrow-1][0..ncol-1] */
float **matrix(long nrow, long ncol)
{
  int i;
  float **m;
  m=malloc_check(nrow*sizeof(float *));
  if ((*m=calloc(nrow*ncol,sizeof(float)))==NULL) {
    fprintf(stderr,"Error allocating %ld by %ld matrix\n",nrow,ncol);
    perror("calloc"); exit(EXIT_FAILURE);
  }
  for(i=1;i<nrow;i++) m[i]=m[i-1]+ncol;
  return m;
}

float ***tensor3(int n1, int n2, int n3)
{
  int i, j;
  float ***m;
  m=(float ***)malloc_check(n1*sizeof(float **));
  *m=(float **)malloc_check(n1*n2*sizeof(float *));
  if ((**m=(float *)calloc(n1*n2*n3,sizeof(float)))==NULL) {
    fprintf(stderr,"Error allocating %d * %d * %d tensor\n",n1,n2,n3);
    perror("calloc"); exit(EXIT_FAILURE);
  }
  for(i=0;i<n1;i++) {
    if (i) {m[i]=m[i-1]+n2; *(m[i]) = *(m[i-1]) + n2*n3;}
    for (j=1; j<n2; j++) m[i][j] = m[i][j-1]+n3;
  }
  return m;
}

float complex ***ctensor3(int n1, int n2, int n3)
{
  int i, j;
  float complex ***m;
  m=(float complex ***)malloc_check(n1*sizeof(float complex **));
  *m=(float complex **)malloc_check(n1*n2*sizeof(float complex *));
  if (!(**m=(float complex *)calloc(n1*n2*(n3/2 + 1),sizeof(float complex)))) {
    fprintf(stderr,"Error allocating %d * %d * %d complex tensor3\n",n1,n2,n3);
    perror("calloc"); exit(EXIT_FAILURE);
  }
  for(i=0;i<n1;i++) {
    if (i) {m[i]=m[i-1]+n2; *(m[i]) = *(m[i-1]) + n2*(n3/2 + 1);}
    for (j=1; j<n2; j++) m[i][j] = m[i][j-1]+(n3/2 + 1);
  }
  return m;
}
// this was commented out but is used in openmp_ic.c (perhaps not used for the two field case anyway); so i am uncommenting it out -saroj (4/15/2013) - but this gives trouble for cross_powerspc2 which has its own fopen_check and I don't need openmp_ic.c so i don't need this (5/16/2013)
/*FILE *fopen_check(char *filename, char *mode)
{
  FILE *fp = fopen(filename,mode);
  if (fp==NULL) {
    fprintf(stderr, "Error opening file %s with mode %s\n", filename, mode);
    perror("fopen"); exit(EXIT_FAILURE);
  }
  return fp;
}*/

size_t ftwrite(void *data, size_t size, size_t nmemb, FILE *fp)
{
  int nbyte = size*nmemb;

  if (nbyte <= 0) {
    fprintf(stderr, "Error in ftwrite: size=%ld, nmemb=%ld, nbyte=%d\n",
	    size, nmemb, nbyte);
    exit(EXIT_FAILURE);
  }
  if (fwrite(&nbyte, 4, 1, fp) != 1) {
    fputs("Error writing to file in ftwrite\n", stderr);
    perror("fwrite"); exit(EXIT_FAILURE);
  }
  if (fwrite(data, size, nmemb, fp) != nmemb) {
    fputs("Error writing to file in ftwrite\n", stderr);
    perror("fwrite"); exit(EXIT_FAILURE);
  }
  if (fwrite(&nbyte, 4, 1, fp) != 1) {
    fputs("Error writing to file in ftwrite\n", stderr);
    perror("fwrite"); exit(EXIT_FAILURE);
  }
  return nmemb;
}



// random complex number where real & imag parts are Gaussian distributed
float complex crangauss()
{
  double y;
  do y = drand48(); while (y<1e-80);
  return sqrt(-log(y))*cexp(2 * M_PI * I *drand48());
}

// Numerical Recipes Gaussian random number generator
float gasdev() 
{
  static int iset=0;
  static float gset;
  float fac,rsq,v1,v2;

  if  (iset == 0) {
    do {
      v1=2.0*drand48()-1.0;
      v2=2.0*drand48()-1.0;
      rsq=v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    fac=sqrt(-2.0*log(rsq)/rsq);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  } else {
    iset=0;
    return gset;
  }
}

int get_option(int argc, char *argv[], char format[],...)
{
  va_list ap;
  char optionstr[100];
  char *tmp;
  int i;
 
  strcpy(optionstr,format); /* optionstr points at the option */
  tmp=strpbrk(optionstr," \t\n");
  if (tmp!=NULL) {
    tmp[0]='\0';
    tmp++;
    tmp=strchr(tmp,'%');              /* tmp points at the format string */
  }
  /*  printf("|%s| |%s|\n",optionstr,tmp); */
  i=0;
  while ((i<argc)&&(strcmp(argv[i],optionstr))) i++;
  if (i==argc) return 0;                  /* option doesn't exist */
  if ((i<argc) && (tmp==NULL)) return 1;  /* option found and not expected */
                                          /* to have arguments */
  /* advance to the first value to be read */
  optionstr[0]='*';                       /* null out the option */
  va_start(ap,format); /* start with first variable address */
 
  ++i;
  /*  printf("|%s| |%s|\n",argv[i],tmp);*/
  while ((i<argc) &&
         (tmp!=NULL) &&
         ((argv[i][0]!='-') || (isdigit(argv[i][1]))))
    {
  
/*       printf("|%s| |%s|\n",argv[i],tmp); */
      switch (tmp[1]) {
 
      case 'd':
        sscanf(argv[i],"%d",(va_arg(ap,int *)));
        break;
      case 'f':
        sscanf(argv[i],"%f",(va_arg(ap,float *)));
        break;
      case 's':
        sscanf(argv[i],"%s",(va_arg(ap,char *)));
        break;
      case 'l':
        if (tmp[2] == 'f') sscanf(argv[i],"%lf",(va_arg(ap,double *)));
        else if (tmp[2] == 'd') sscanf(argv[i],"%ld",(va_arg(ap,long *)));
        else return -1;
        break;
      default: return -1; break;
      }
 
      i++;
      tmp++;
      tmp=strchr(tmp,'%');
    }
 
  va_end(ap);
 
  if (i==argc) return -1;
 
  return 1;
}
